home *** CD-ROM | disk | FTP | other *** search
/ AP Professional Graphics CD-ROM Library / AP Professional Graphics CD-ROM Library.iso / pc / unix / appendix / gemsiii / urot.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-03-13  |  3.9 KB  |  116 lines
  1. /* urot.c */
  2. /* Generates a uniform random rotation */
  3. /* Ken Shoemake, September 1991 */
  4.  
  5. #include <stdlib.h>
  6. #include <math.h>
  7. #include "GraphicsGems.h"
  8.  
  9. /* Define an INT32 value to be a 32 bit signed integer */
  10. typedef int INT32;
  11.  
  12. typedef struct {float x, y, z, w;} Quat;
  13. enum QuatPart {X, Y, Z, W, QuatLen, V=0};
  14.  
  15. /* * * * * *  Utility for quaternion conversion * * * * * */
  16.  
  17. /** Qt_ToMatrix
  18.  *  Construct rotation matrix from quaternion (unit or not).
  19.  *  Assumes matrix is used to multiply row vector on the right:
  20.  *  vnew = vold mat.  Works correctly for right-handed coordinate system
  21.  *  and right-handed rotations. For column vectors or for left-handed
  22.  *  coordinate systems, transpose the matrix.
  23.  */
  24. void Qt_ToMatrix(Quat q, Matrix3 *out)
  25. {
  26.     double norm = q.x*q.x + q.y*q.y + q.z*q.z + q.w*q.w;
  27.     double s = (norm > 0.0) ? 2.0/norm : 0.0;
  28.     double xs = q.x*s,        ys = q.y*s,        zs = q.z*s;
  29.     double wx = q.w*xs,    wy = q.w*ys,    wz = q.w*zs,
  30.              xx = q.x*xs,    xy = q.x*ys,    xz = q.x*zs,
  31.              yy = q.y*ys,    yz = q.y*zs,    zz = q.z*zs;
  32.     double (*mat)[3] = out->element;
  33.     mat[X][X] = 1.0 - (yy + zz); mat[X][Y] = xy + wz; mat[X][Z] = xz - wy;
  34.     mat[Y][X] = xy - wz; mat[Y][Y] = 1.0 - (xx + zz); mat[Y][Z] = yz + wx;
  35.     mat[Z][X] = xz + wy; mat[Z][Y] = yz - wx; mat[Z][Z] = 1.0 - (xx + yy);
  36. } /* Qt_ToMatrix */
  37.  
  38.  
  39. /* * * * * *  How to do it using gaussians * * * * * */
  40.  
  41. /** Qt_RandomG
  42.  *  Generate uniform random unit quaternion from random seed.
  43.  */
  44. Quat Qt_RandomG(INT32 *argseed)
  45. {
  46. /*  This algorithm generates a gaussian deviate for each coordinate, so
  47.  *  the total effect is to generate a symmetric 4-D gaussian distribution,
  48.  *  by separability. Projecting onto the surface of the hypersphere gives
  49.  *  a uniform distribution.
  50.  */
  51.     Quat q;
  52.     /* uurand generates doubles uniformly distributed between ╨1 and +1 */
  53.     /* This linear congruential generator is inline to exploit signed ints */
  54.     register INT32 seed = *argseed;
  55. #define uurand() ((seed = (seed+1)*69069)/2147483648.0)
  56.     register double x = uurand(), y = uurand();
  57.     register double z = uurand(), w = uurand();
  58.     register double s1, s2;
  59.     double num1, num2, root1, root2, r;
  60.     while ((s1 = x*x+y*y) > 1.0) {x = uurand(); y = uurand();}
  61.     while ((s2 = z*z+w*w) > 1.0) {z = uurand(); w = uurand();}
  62.     /* Now the point (x,y) is uniformly distributed in the unit disk */
  63.     /* So is the point (z,w), independently */
  64.     num1 = -2*log(s1); num2 = -2*log(s2);
  65.     /* Now x*sqrt(num1/s1) is gaussian distributed, using polar method */
  66.     /* Similarly for y, z, and w, and all are independent */
  67.     r = num1 + num2;    /* Sum of squares of four gaussians */
  68.     root1 = sqrt((num1/s1)/r); root2 = sqrt((num2/s2)/r);
  69.     /* Normalizing onto unit sphere gives uniform unit quaternion */
  70.     q.x = x*root1; q.y = y*root1; q.z = z*root2; q.w = w*root2;
  71.     *argseed = seed;
  72. #undef uurand
  73.     return (q);
  74. } /* Qt_RandomG */
  75.  
  76. /** M3_RandomRotG
  77.  *  Generate uniform random rotation matrix from random seed.
  78.  */
  79. void M3_RandomRotG(INT32 *seed, Matrix3 *m)
  80. {
  81.     Qt_ToMatrix(Qt_RandomG(seed), m);
  82. } /* M3_RandomRotG */
  83.  
  84.  
  85. /* * * * * *  How to do it using subgroup algorithm * * * * * */
  86.  
  87. /** Qt_Random
  88.  *  Generate uniform random unit quaternion from uniform deviates.
  89.  *  Each x[i] should vary between 0 and 1.
  90.  */
  91. Quat Qt_Random(double x[3])
  92. {
  93. /*  The subgroup algorithm can be condensed to this efficient form.
  94.  *  Use rotations around z as a subgroup, with coset representatives
  95.  *  the rotations pointing the z axis in different directions.
  96.  */
  97.     Quat q;
  98.     register double r1 = sqrt(1.0 - x[0]), r2 = sqrt(x[0]);
  99.     register double t1 = PITIMES2*x[1], t2 = PITIMES2*x[2];
  100.     register double c1 = cos(t1), s1 = sin(t1);
  101.     register double c2 = cos(t2), s2 = sin(t2);
  102.     q.x = s1*r1; q.y = c1*r1; q.z = s2*r2; q.w = c2*r2;
  103.     return (q);
  104. } /* Qt_Random */
  105.  
  106. /** M3_RandomRot
  107.  *  Generate uniform random rotation matrix from uniform deviates.
  108.  */
  109. void M3_RandomRot(double x[3], Matrix3 *m)
  110. {
  111.     Qt_ToMatrix(Qt_Random(x), m);
  112. } /* M3_RandomRot */
  113.  
  114. /* End of urot.c */
  115.  
  116.